The Correlation Energy and the Spin Susceptibihty of the Two- Valley 

Two-dimensional Electron Gas. 
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We find that the spin susceptibility of a two-dimensional electron system with valley degeneracy 
does not grow critically at low densities, at variance with experimental results [A. Shashkin et at, 
Phys. Rev. Lett. 96, 036403 (2006)]. We ascribe this apparent discrepancy to the weak disorder 
present in experimental samples. Our prediction is obtained from accurate correlation energies 
computed with state of-the-art diffusion Monte Carlo simulations and fitted with an analytical 
expression which also provides a local spin density functional for the system under investigation. 
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I. INTRODUCTION 

The spin properties of low-dimensional electron sys- 
tems in solid state devices are of great interest in re- 
lation to spintronics and quantum computing^, both 
at the fundamental level and for technological applica- 
tions, the long wavelength spin susceptibility of the two- 
dimensional electron gas (2DEG) playing an important 
role in the control of nuclear spins^. They are also be- 
lieved to be intimately related to the apparent metal- 
insulator transition (MIT) observed in 2D'^''*'^'^'^. Indeed 
the spin susceptibility Xs of the 2DEG, measured with 
various techniques^, is consistently found to grow with re- 
spect to its noninteracting Pauli value xo , as the density 
is lowered and the MIT approached^iii^. Recently, exper- 
imental evidence has been given for a critical growth of 
Xs in Si-MOSFETs at a finite density^ coincident, within 
experimental uncertainties, with the critical density for 
the MIT^'^. The qualitative question to which we give 
an answer in this paper is whether such a divergence is a 
property of the ideally clean two- valley (2V) 2DEG, the 
simplest model of electrons confined in a Si-MOSFET— , 
or is due to some other factor. It should be stressed 
from the outset that the valley degree of freedom has 
qualitative effects on the 2DEG properties, making the 
fully spin polarized fluid never stable^, at variance with 
the one-valley (IV) 2DEG, and importantly affects the 
MITiSiii. 

Correlation plays a crucial role in the so-called EG, i.e., 
electrons with a 1/r pair potential, moving in a neutral- 
izing charge backgrounds^. Its importance grows both 
with lowering the density and the space dimensional- 
ity and tends to quantitatively and often even qualita- 
tively change the predictions of simple schemes, such as 
the Hartree-Fock (HF) or the random-phase approxima- 
tion (RPA)i^. In the low-density strongly-correlated EG, 
which would be more properly called an electron liquid, 
the energy balance determining the system properties is 
played on a very minute scale and, to get meaningful 
predictions, a great accuracy such as the one afforded by 
quantum Monte Carlo (QMC) methods is necessaryi^. 



QMC simulations have provided over the years 
the method of choice for microscopic studies of the 
2DEG^iHiiia^ii£a2, which recently has been shown to 
provide a rather accurate model for electrons confined 
in solid state devicesi^. However, no QMC prediction 
is available for Xs in the 2V2DEG and other theoretical 
estimates, obtained either in RPA^^i?? or with a classical 
mappingSl, do not appear reliable^. Here, to calculate 
Xs we resort to extensive state-of-the-art simulations of 
the 2V2DEG, using the diffusion Monte Carlo (DMC) 
technique^^. We thus obtain for the first time the de- 
pendence of the ground state energy on both the density 
and the spin polarization, also improving on Ref. with 
the use of twist-averaged boundary conditions (TABC)^ 
and trial wavefunctions including backflow (BF}i^. 



II. CORRELATION ENERGY OF THE 2V2DEG 

In the 2V2DEG electrons possess an additional discrete 
degree of freedom, i.e. the valley flavor or index, which 
can be conveniently described with a pseudospin. One 
may identify electrons with given spin and pseudospin 
indexes as belonging to a species or component. Accord- 
ingly, the paramagnetic 2V2DEG is a four-component 
system, while both the fully spin-polarized 2V2DEG and 
the paramagnetic 1V2DEG have two components. For 
the sake of simplicity, we restrict here to the symmetric 
case where the number of electrons and the spin polariza- 
tion are the same for both valleys^*^. Thus, at zero tem- 
perature, the state of the system is fully specified by the 
coupling parameter Vs = \/ ^-k naB and the spin polar- 
ization C, = (n|— n|)/n, with n the total electron density, 
Ob the Bohr radius, n-^(i) the density of up (down) spin 
electrons. Below, Rydberg units are used throughout. 



A. Simulation details 

We have performed simulations with the fixed-phase 
(FP)^ DMC method, which gives the lowest upper 
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TABLE I; Parameters of the analytic representation ((2)1 and 
([3]) of the correlation energy of the 2V2DEG, determined from 
Eq. IT} by a least squares fit to the data listed in Table IIIII 
The reduced chi square is = 4.82. A * marks constrained 
parameters, whereas C2 is fixed to zero since it turned out to 
be irrelevant in the fitting procedure. The parameters rj, rj^, 
7 and 72 in Eq. ((T} only concern the size extrapolation; their 
optimal values are 0.056, 0.17, 2.03 and 0.45, respectively. 





i = 


i = 1 


i = 2 


A, 


-0.99870* 


0.44570* 


0.0082290 




if(10-3^)* 


-0.85288* 


0.048979 




0.62208 


-7.6202 
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-1.6194 


-0.051302 


E, 


1.6208 


12.714 


25.911 


F, 


-0.012856 


0* 


0* 


G, 


0.66150 


19.692* 


15.072* 




0.029765* 


3.6334* 


6.2343* 


P 




11.879 





bound to the ground-state energy consistent with the 
many-body phase of a suitably chosen, complex-valued 
trial function. For real trial functions FP-DMC reduces 
to the standard fixed-node (FN) approximation^^. A 
complex trial function allows using TABC^'i, which re- 
duce the size dependence of the kinetic energy by one 
order of magnitude with respect to periodic boundary 
conditions (PBC). Furthermore, since TABC do not re- 
quire closed shells in /c-space there are no restrictions 
on the number of electrons per component, so that the 
polarization can be changed by flipping any number of 
spins, with fixed total number of electrons^l. Our trial 
function is the product of Slater plane-wave (PW) de- 
terminants (one per component) and a Jastrow factor—. 
BF correlations^ are included only for ^ = and C = 1; 
but with FN-DMC and in PBC. Their contribution to 
the ground-state energy is then added to the PW en- 
ergies assuming a quadratic dependence on polarization 
as m Ref . miia. The ground-state energy per particle 
EN{rs,C): calculated for several values of Vg, (, and the 
electron number TV, is recorded in Table Hill of Appendix 
A. 



B. Analytic representation 

Following Ref. we determine the energy per particle 
E{rs, C) the thermodynamic limit by fitting to the data 
listed in Table Hill an analytic expression which embodies 
the Vs and ^ dependence as well as a Fermi-liquid-like 
size correction: 
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+ (7 + 7.0 7- 



TsN 

(1) 



The fitting parameters rj,riz,^ and 72 take into account 
potential energy finite-size effects, while AT/v(r5, C) is the 
difference of the non-interacting kinetic energy evaluated 
at finite N with TABC and in the thermodynamic limit. 
E{rs,^) is customarily decomposed as sum of the non- 
interacting kinetic energy, eo(?'s, C) = (1 + C^)/(2rs), the 
exchange energy, e,(r„C) = -^;^[(l+C)^/'+(l-C)^/'], 
and the unknown correlation energy ec(rs,C), for which 
we adopt the same analytical representation of Ref. [TtI . 



l)ei«)(r„C)+ V C"«.(r-.), (2) 
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where ei''(r„ C) = e,(r„ C) - (1 + |C' + ^C^)eArs,0), 
and the functions ai{rs) are defined by 

1 \ 
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(3) 



We constrain the correlation energy ^ to satisfy known 
high- and low-density limits (Appendix B), reducing in 
this way the number of free fitting parameters from 29 to 
18. The correlation energy of the 2V2DEG, as given by 
Eqs. ([2]) and ([3]) with the parameters listed in Table [H 
represents a central result of this work. 



C. Phase diagram 

In Fig. [1] we plot the energies of the paramagnetic 
and the fully spin-polarized 2V2DEG. They are shown 
by solid lines labeled with the number of components, 
Nc — 4: and Nc = 2, respectively. For comparison, we 
also plot QMC results for other phases of the 2DEG: 
A'c = 1 labels the fully polarized one- valley 2DEG^-, 
whereas the 2D charged-boson fluid^^ corresponds to the 
limit of an infinite-component 2DEG. The energy of the 
Wigner crystal is known to be almost independent of 
the number of component a^'^i^^ ; we report here the re- 
sult of Ref. [1^. We note that at large as the number 
of components increases the ground-state energy appears 
to quickly approach the infinite-component limit. 

The dashed line for iVc = 2 in Fig. [Tj is the result of 
Ref. [l3 for the paramagnetic 2DEG. Its agreement with 
our curve for the polarized 2V2DEG is expected, but still 
gratifying: the two calculations differ by details in the ex- 
trapolation to the thermodynamic limit, and the close- 
ness of their results supports a good control of the finite- 
size bias. The dashed line for iVc = 4 is instead the result 
of Ref. 1^ for the paramagnetic 2V2DEG. Its difference 
with the present results conveys a physical information, 
namely the quantitative effect of BF correlations which 
were not included in the previous simulations^ . Backflow 
improves the nodal structure of the PW wave function, 
yielding in the FN approximation a tighter upper bound 
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FIG. 1: Energy per particle of various phases of the 2DEG at 
T = 0. The energy label indicates the number Nc of equiv- 
alent components in the homogeneous fluid, or the bosonic 
fluid (B), or the triangular Wigner crystal (WC), as applica- 
ble. Results of the present study are given by solid lines. The 
dashed lines are from Ref.[i3 ( Ei, E2 ), Ref.fi (E4 ), Ref.[ll 
(Ewe), and Ref.IH (Eb ). The inset shows E{rs, C)-E{rs,0) 
from Eq. ([2]) (solid line) together with the simple quadratic 
dependence [E{r,, 1) - E{rs, 0)]C (dotted), for the 2V2DEG, 
and the result for the 1V2DEG^ (dashed) magnified by a 
factor 10, at — 25. 



to the exact ground-state energyA^. It is knowiJ^ that 
BF correlations lower the FN energy more for Nc = 2 
than for iVc = 1. Here we find that the BF energy gain 
for Nc = 422i is smaller than for Nc — 2 (see Table HI] 
of Appendix A), albeit larger than for Nc = 1. The 
modest effect of BF correlations for Nc — 4 entails only 
marginal quantitative changes to the phase diagram of 
the 2V2DEG predicted in Ref. i. The density of Wigner 
crystallization shifts to a slightly lower value, ~ 45. 

Before discussing the spin polarization dependence of 
the energy and our prediction for the spin susceptibility, 
we should stress that our results provide the most ac- 
curate available estimate for the correlation energy Cc of 
the 2V2DEG, which in turn is the key ingredient for den- 
sity functional theory (DFT) studies of inhomogeneous 
two-valley systems in 2D within the local spin density 
approximation^. The knowledge of Cc allows also to 
check the accuracy of the ansatz made in Ref. to con- 
struct the correlation energy of a system with an arbi- 
trary number of components, edNc), in terms of that of 
the one valley systenJ^. A comparison between edNc) 
from Ref. i31. . the present Cc, and the nominally exact 
QMC results for charged bosonsS^ (Appendix C) exposes 
the limited accuracy of edNc) especially at large r^, in- 
cluding its prediction— of an unphysical transition be- 
tween the nodeless ground state of the infinite-component 
system and the antisymmetric ground state of the one- 



FIG. 2; Spin-susceptibility enhancement of the 1V2DEG and 
the 2V2DEG. The results of the present work are compared 
with HF and RPAi2ii.o pj-edictions as well as QMC results for 
the one valley casei^. Experimental results for Si-MOSFETs 
are also shown^i^. 



component 2DEG. Yet, the comparison between DFT 
calculations of two- valley symmetric systems using either 
^d-^c) or the present Cc would provide a valuable test of 
the adequacy of ec{Nc) for DFT applications. 

Our calculations confirm the absence of a transition 
from the paramagnetic to the fully spin-polarized fluid 
in the two-valley symmetric system^. Moreover, in the 
whole density range where the fluid is stable we find no 
evidence for the stability of a state with partial spin po- 
larization. As illustrated in Fig. [T] for = 25, E{rs,() 
displays its minimum at C = and, for all practical pur- 
poses, can be considered a convex function of Con- 
vexity ensures that, by turning on an in- plane magnetic 
field B, the absolute minimum displayed by the energy 
goes continuously from C = to C = 1- If the energy 
exhibits a local maximum or even an inflection point for 

< 1, instead, the B-driven transition to the full spin 
polarization becomes a first-order one and is accompa- 
nied by a jump in the polarization^. This is clearly the 
case for the 1V2DEG at rg = 25 also shown in Fig. [TJ 



III. SPIN SUSCEPTIBILITY 

The spin susceptibility enhancemenl^ of the 2V2DEG 
is readily calculated using Eq. ^ as 



Xs/xo 



1 



2r., 



2r2ai(r,) 



(4) 



In Fig. [2] we compare our QMC prediction with the avail- 
able experimental results for electrons confined in Si- 
MOSFETs. It is evident that the 2V2DEG spin sus- 
ceptibility moderately overestimates experiments at high 
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density but largely underestimates them at low density, 
where it does not display any critical growth. In fact Xs 
is a concave function of Vg at all densities where the fluid 
phase is stable. Indeed, a realistic description of a 2DEG 
in a solid state device requires consideration of additional 
elements such as transverse thicknes o^^'^^ and disorder 
scattering"-. As the thickness is known to suppress the 
spin susceptibility and a weak disorder to enhance it, at 
present the only likely candidate to explain the experi- 
mentally observed critical behavior of appears to be a 
weak disorder. In Fig. [5] we also report the QMC results 
of a 1V2DEG^^. It is clear that the valley degeneracy 
causes a substantial suppression of the spin susceptibility, 
in qualitative agreement with the effect found in experi- 
ments on AlAs based quantum wells^, though for an in- 
plane anisotropic mass. Moreover, Xs(^s) changes from 
a convex to a concave function in going from the IV to 
the 2V2DEG. We also show in the figure the predictions 
of HF and RPA. The general trend is that, while RPA 
performs somewhat better than HF, both largely overes- 
timate the QMC predictions and yield divergences which 
either have no counterpart in QMC, for the 2V2DEG, or 
in the best case take place at a density about 13 times 
larger than in QMC, for the 1V2DEG. At least RPA 
reverts the qualitatively wrong prediction of HF which 
yields an enhancement of the spin susceptibility in going 
from the one- to the two-valley system. 



IV. DISCUSSION AND CONCLUSIONS 

We have reliably estimated the spin susceptibility of 
the 2V2DEG, which provides the simplest model for elec- 
trons confined in Si-MOSFETs. Our results clearly point 
to the crucial, qualitative role of weak disorder scattering 
in determining the critical growth found in the measured 
susceptibility at low density^ and to a likely minor, quan- 
titative role of transverse thickness in suppressing the 
susceptibility at high density. 2D electron systems in high 
mobility Si-MOSFET's at times have been termed clean, 
meaning in fact without admixture of local moment^, but 
also implicitly implying that observed properties would 
be disorder independent and would correspond to those 
of an ideally clean electron gas. This latter viewpoint, 
fostered by the recent experimental observation that the 
effective mass enhancements of samples^'''^^ with peak 
mobilities differing by about one order of magnitude ap- 
pear to be the same within error bars (of about 10%), 
is contradicted by our findings. We should stress indeed 
that the samples of Refs. [SJ and [s^ are different on a 
number of counts and not only for the amount of disor- 
der. Electrons in (111) Si-MOSFETs^ have (i) a sizeable 
band mass anisotropy nix/rny = 0.28, i.e. comparable 
with the one in AlAs quantum wells^, and (ii) a trans- 
verse thickness parameter v^3)6/(rsa_B) (see, e.g., Ref.fisl 
for the definition) which is more than twice the one in 
(100) Si-MOSFETs^^. Both effects (mass anisotropy and 
thickness) are known to suppress spin susceptibility in an 



appreciable manner-isi^. Moreover, comparing the abso- 
lute peak mobilities of Ref. [H on the one hand and of 
Ref. |35j on the other, i.e. of systems with quantitatively 
different length and energy scales (due to different band 
masses) is not appropriate. If I and as are respectively 
the mean-free path and the effective Bohr radius in a 
given system, we find that the peak of l/{rsaB) for the 
EG of Ref. i34 is only 3 times smaller than that of the 
EG of Ref. [3^. Hence the experiment in Ref. [H in our 
opinion is not at all conclusive in ruling out an effect 
of disorder on the effective mass, let alone on the spin 
susceptibility of these systems. 

We have also obtained: an assessment of the back- 
flow effects on the energy of the two- valley paramagnetic 
phase, which remains stable with respect to any partially 
or fully polarized phase, up to the Wigner crystallization; 
an analytical fit of the QMC correlation energy, which 
also interpolates between exact high and low-density lim- 
its, and provides a local spin density functional for DFT 
studies of two-valley systems; the clear indication that 
an accurate account of correlation beyond RPA is cru- 
cial when considering the properties of both the IV- and 
2V2DEG. 



APPENDIX A: DETAILS OF THE DMC 
SIMULATIONS 

The trial function was chosen of the usual Slater- 
Jastrow form, "^{K) = D{R)J{R), where R = 
(ri,...,rAr) represents the coordinates of the N elec- 
trons. The Jastrow factor is a pair product, J{R) = 
exp [ — X]i<j ""(^u )] : with u{r) the parameter-free RPA 
pseudopotentiali^. The phase structure is fixed by the 
complex factor D = Yiv^t^^ i-^-' ^ product of Slater de- 
terminants, one for each spin- valley component. 

Most of the simulations were carried out with the stan- 
dard plane-wave (PW) choice for the one-particle or- 
bitals, = det[exp(ik, • r^)]. For C = and C = 1 

we also included backflow (BE) correlational^, D^^ = 
det[exp(ik^ • x^)], where = -I- J2j^iV{nj){'r^ - Vj) 
and the BF function ri{r) (of the form suggested ir^ii) 
was optimized by minimization of the variational energy. 

We simulated the imaginary-time evolution of the sys- 
tem by a branching random walk, using a short-time ap- 
proximation of the importance-sampled Green's function 
and exerting control on the number of walkers. Calcu- 
lations were performed at Vs — 1, 2, 5, 10, 20, 40. For 
C = and C = 1 we chose several values of the number 
of electrons between = 36 and N = 116, whereas 
11 intermediate values of the polarization, defined by 
flipping one spin at a time, were studied for N = 52. 
The twist average, for the PW simulations, was per- 
formed on a mesh defined by qx{i) = A(* — 1/2), QyU) = 
A{j - 1/2), 1 < i < 8, i < j < 8, A = tt/SL, with L the 
side of the simulation box. Long-range interactions were 
dealt with the optimized-splitting method of Ref. [stI 

Extrapolation to zero time step r and infinite number 
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TABLE II: Difference A = E§''{r,,Q - E^^{rs,0 between 
the BF and the PW energy (in Rydberg per particle) at se- 
lected values of Va, A'^. In parentheses the statistical error 
on the last digit. 





C = o 


C = i 




TV A 


iV A 


1 


52 -0.0028(1) 


50 -0.0034(1) 
58 -0.0035(2) 
90 -0.0032(1) 


2 


52 -0.00166(5) 


42 -0.00175(9) 
50 -0.00192(9) 
58 -0.00217(9) 


5 


52 -0.00057(2) 


50 -0.00077(3) 
58 -0.00088(3) 


10 


52 -0.00021(1) 
84 -0.00022(1) 


42 -0.00025(1) 
50 -0.00030(2) 
58 -0.00032(2) 


20 


52 -0.000043(6) 


42 -0.000081(7) 
50 -0.000085(7) 
58 -0.000116(6) 
90 -0.000116(6) 


40 


52 -0.000020(1) 


50 -0.000020(1) 
90 -0.000031(1) 



of walkers Nw was also carried out at fixed density, on 
the assumption that the r and Nw dependences are ap- 
proximately independent. Results at polarizations C = 0; 
( ~ 0.5, and C = 1 and for a bunch of r (Nw) values were 
used to establish the r ( Nw) dependence of the energy 
as function of these dependences, combined together, 
were then used to extrapolate to Nw = oo, r = the 
energies calculated for all values of 

We record the difference between BF and PW energies 
at zero and full polarization in Table |TT] and the whole 
set of energies extrapolated to Nw = cxd, t = and 
including the backflow correction in Table Hill 

APPENDIX B: HIGH AND LOW DENSITY 
LIMIT OF THE CORRELATION ENERGY OF 
THE FOUR COMPONENT 2DEG 

We directly refer to Ref. [l^for both the ^ oo limit, 

^ 3/2 

whose leading terms in r~ and are independent 

of the number of components^, and the ^ limit at 
^ = 1, which is the same two-component system as the 
one- valley case at ( = 0^. 

Here we only need to specify the high-density limit 
for the four-component system, lim^.^p ec(rs,0) = Aq + 
Bq rg In r s . Generalizing the procedure of Ref . |39| to the 
multivalley case, we write ec as the sum of the second- 
order exchange energy ej'''' and the ring contribution ei^"^ , 



whose lowest order is the direct term of the second- 
order energy per particle. It turns out that £3^^ is a con- 
stant, independent of and the number of components 




10 20 30 40 

FIG. 3: Phase diagram of the multicomponent 2DEG: liquid 
phases. Solid E4 and E2 from the present work, solid Ei and 
long-dashed E2 from Ref. [13, Eb from Ref. HI. £4* and Eb- 
are from Ref. |3l|. One- and two-component energies from 
Ref. [S^ coincide with Ref. ones by construction. 



of the system, Nc, while = (-^c) depends only on 
Nc- Furthermore, we notice that the ring contribution 

scales with Nc as ei^\rs, Nc) = Nc f{rs Nc^'^), so that 
the following scaling law holds 

eW(r„4) = 2eW(4\/2r„2). (Bl) 

By applying the scaling law (|B1[) to the leading terms 

of ec, we find Aq = e^''^ -I- 2e^'^^(2) = -0.99870, Ba = 
16(10- 37r)/(37r). 



APPENDIX C: CHECK OF AN APPROXIMATE 
MULTICOMPONENT CORRELATION ENERGY 

In Fig. [3] we show a comparison between the multicom- 
ponent correlation energy ec{Nc) of Ref. [3l| and various 
simulation results, including the present two- valley calcu- 
lation, and the nominally exact QMC results for charged 
bosons^*. Total energies are displayed. The scale of the 
figure, emphasizes the limited accuracy of ec{Nc) in the 
large regime. 



^ Semiconductor Spintronics and Quantum Computing, edited by D. D. Awscholom, N. Samarth, and D. Loss 
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TABLE III: Data used for the fit described in tlie paper. Twist-averaged DMC energy in Rydberg per particle i5jv(r-a, (^), 
calculated at finite A*', extrapolated to zero time step and infinite number of walkers, and including BF correlations; in 
parentheses the statistical error on the last two figures shown. The backflow correction was obtained from Table HIl employing 
the results at the largest A'' available. 
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52 
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52 


11/13 


-0.301915(42) 




52 


10/13 
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52 
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52 
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52 
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84 
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84 
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